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1 Discussion 



in 

o. 

Dr Sakov cited an earlier work |llj . which is similar to [6], but whose con- 
clusion contradicts that in j6]. We note that these conclusions shall be inter- 
preted based on how the filter is implemented. In [11] , the authors considered 
the combination of the ensemble transform Kalman filter (ETKF) with the 
positive-negative pair (PNP) scheme. In order to keep the size of the analysis 
ensemble constant at each assimilation cycle, the authors in [H] choose to 
discard half of the members before propagating them forward (otherwise the 
ensemble size will double at each cycle). The advantage in doing this is that 
the computational cost can be very cheap. However, as a potential problem, 
the sample covariance of the remaining ensemble members might not be a good 
approximation to the one before halving the size, which may thus deteriorate 
the performance of the filter. 
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Dr Sakov also compared the numerical results of [6] with the existing ones in 
the literature, for example, [Tf9][T2] . In terms of absolute rms error, Dr Sakov 
commented that the results in jB] are in the range of 0.73 — 0.99 [5], while those 
in [7TI9"][T2"] can be as low as 0.17 — 0.21. This difference may depend on how the 
comparison is made. For example, in order to obtain a better performance, in 
[7] the authors used a bank of ETKFs, while in [9] the ensemble size was up to 
91 for the 40-dimensional Lorenz-Emanuel model (L40 after j8]). In contrast, 
in [6] the number of sigma point is at most 13, since the upper bound l u = 6. 
If in [9] one chose the ensemble size to be 13 for the L40 system, then it can 
be seen that the rms errors of the various ETKFs are all above 1 (cf. Fig. 3 
of [9]), higher than the range of 0.73 — 0.99 achieved by the EnUKF in [6]. 

In [6], we have focused on the scenario where the ensemble size (or the rank 
of the sample covariance) in the filter is typically (much) lower than the di- 
mension of the system in assimilation. Nevertheless, it should be interesting to 
compare the performances of the EnUKF and the EnKF (with the ETKF as 
the representative) beyond the aforementioned scenario, where the ensemble 
size of a filter could be close to, or even higher than, the dimension of the 
system in assimilation. To this end, in what follows we would like to conduct 
one more experiment, supplementary to those in [5]. 

In the experiment, we examine the performances of the EnUKF and the ETKF 
under various covariance factors (but without covariance filtering to avoid 
complicating our discussion). The settings in the experiment are as follows. 
The truths and observations are generated in the same way as that in |JJ 
§ 4]. For the EnUKF, we let its intrinsic parameters (3 = 2, A = —2, the 
threshold h\ = 1000. We let the lower bound l\ equal to the upper bound 
l u (hence the truncation number I — h — l u ), and take their values in the 
set [5,10,15,20,30,40], which implies that the corresponding ensemble size 
n — 21 + 1 takes values in the set [11, 21, 31, 41, 61, 81]. In addition, we let the 
inflation factor S take values in the set [[0.05 : 0.05 : 0.5], [1 : 0.5 : 6]]. Here, the 
notation [0.05 : 0.05 : 0.5] denotes the set whose elements increase from 0.05 
to 0.5, with an even increment of 0.05 each time. Similar interpretation can 
also be applied to the notation [1 : 0.5 : 6], while [[0.05 : 0.05 : 0.5], [1 : 0.5 : 6]] 
means the concatenation of the sets [0.05 : 0.05 : 0.5] and [1 : 0.5 : 6]. For 
comparison, we use the same ensemble sizes n and inflation factors 5 in the 
ETKF. At the first assimilation cycle, the EnUKF is always initialized with a 
10-member background ensemble, no matter what the value of I (hence n) is. 
But for the ETKF, it is initialized with an n-member background ensemble. 
To reduce the statistical fluctuation, we repeat the experiment for both the 
EnUKF and the ETKF for 20 times, each time with the same initial system 
state (hence the subsequent truths) and observations, but a randomly drawn 
initial background ensemble. 

In Figs. [T]and [2J we show the relative rms errors of the EnUKF and the ETKF 
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Fig. 1. Relative rms errors of the EnUKF (with different ensemble sizes) as functions 
of the inflation factors. 




Fig. 2. Relative rms errors of the ETKF (with different ensemble sizes) as functions 
of the inflation factors. 

(with different ensemble sizes) as functions of the inflation factor. As one can 
see in them, in order to achieve a relatively good performance, the choice of 
the inflation factor may be different with different ensemble sizes. In the large 
sample scenario (say, n > 40), a small covariance inflation factor, e.g., 5 = 0.05, 
appears good for both the EnUKF and the ETKF, which is consistent with 
Dr Sakov's comment. However, in the small sample scenario (say, n = 11), 
one may need a larger inflation factor, e.g., 5 = 3 for the ETKF and 5 = 6 for 
the EnUKF, to achieve a relatively good result. Our explanation is that in the 
small sample scenario, there exist relatively large sampling errors. Hence, a 
larger inflation factor is needed to compensate for the effect of finite ensemble 
size. 



Given the same ensemble size n, in terms of the lowest relative rms errors that 
the filters can achieve within the tested parameter ranges, the performances 
of the EnUKF and the ETKF are comparable when n = 11, with the lowest 
rmse of the EnUKF being 0.490, slightly lower than that of the ETKF (which 
is 0.493). As n increases to n — 21, the ETKF becomes to outperform the 
EnUKF. If n increases further, the ETKF still outperforms the EnUKF, but 
the gap in their performances are narrowed. For example, the lowest relative 
rmse of the ETKF is 0.049, which is obtained at 5 = 0.05 with n = 61, while 
the lowest relative rmse of the EnUKF is 0.075, achieved at 5 = 0.05 with 
n = 41. 

The above result reflects the difference between the EnUKF and the ETKF 
in their philosophies of statistics estimations. As shown in the Appendix of 
P|, the unscented transform (UT) used in the EnUKF aims to produce the 
estimations which capture the first few terms in the Taylor expansions of the 
true mean of covariance, but which in general mismatch higher order terms. 
In contract, the Monte Carlo approximation adopted in the ETKF produces 
the estimations that converge asymptotically to the true statistics. However, 
given only a small ensemble size, the resulting large sampling errors may lead 
to large biases and spurious modes. For this reason, the relative superiority 
between the EnUKF and the ETKF may invert in the small and large sample 
scenarios. 



2 Conclusion 



Through the supplementary experiment, we showed that, in the small sample 
scenario, the EnUKF can perform (slightly) better than the ETKF with the 
same ensemble size, while as the ensemble size increases, the ETKF outper- 
forms the EnUKF instead. Our explanation of this phenomenon is the follow- 
ing: The estimations of the ETKF are based on the Monte-Carlo approxima- 
tions, which converge asymptotically to the true values as the ensemble size 
increases. However, in the small sample scenario, the Monte-Carlo approxima- 
tions may exhibit large biases and spurious modes due to the relatively large 
sampling errors, which are partially avoided in the estimations based on the 
UT. However, one possible problem of the UT is that its estimations may not 
converge to the true statistics as the ensemble size increases, since its ensem- 
ble members are deterministically chosen, rather than randomly drawn as the 
samples in the Monte-Carlo approximations. 



3 Appendix: Implementation of the covariance filtering technique 

A comparison of the results presented in Figs. CD and [2] (when n = 11) in 
this reply and those in Figs. 2 and 3 of j6j confirms the benefit of covariance 
filtering. By conducting covariance filtering, one increases the rank of a sample 
covariance matrix. So in effect, it increases the effective ensemble size. Other 
benefits (and shortcomings) of covariance filtering are also discussed in, for 
example, [3"fTU] . 

Covariance filtering essentially involves computing the Schur product between 
a sample covariance and a taper matrix. Mathematically, the Schur product, 
C = A o B, of two matrices A and B with the same dimensions, is defined 
as the matrix with the same dimension as A and B, whose components Cij = 
AijBij, where Aij, Bij and CV,- are the components on the ith row and the 
jth column of the matrices A, B and C, respectively. In our discussion, we 
suppose that A is the covariance matrix, and B is the taper matrix used to 
reduce the spuriously large correlations in A. 

The construction of B can be done in the following way: 

By = p(dij), (1) 

where d^ is a metric measuring the difference between the locations % and j, 
and p is a correlation function. Several examples of correlation functions were 
discussed in [2]. In our experiments, we followed [I] and chose function p in 
the following form: 



( 1 S 


1 „ 


5 


5 „ 




— z 5 


+ ^ 


+ -z 3 - 


--z 2 


+ 1 


l\ 


l\ 


$ 3 


§ 2 




—z 5 


--2 4 


+ -z 3 + 


-z 2 


-52 


12 


2 


8 


3 




o, 

V 











if < z < 1 ; 



p(z) = {± z 5_i^ + 2.^ + ^ z 2_ 5 z + 4-± z -i if 1< 2 < 2 ; ( 2 ) 

3 ~ ' 



2 

if z > 2 . 



For illustration, the shape of the function p is plotted in Fig. ()3]). As one can 
see there, the function p has a "cut-off" effect at z = 2, in the sense that the 
values of the function are set to zero for all z > 2. 

For data assimilation in real world, dij is normally a function of the distance 
between the locations % and j in the three dimensional physical world. For 
example, see [T|. However, in mathematical analysis, this choice might not 
always be available. For example, the system states of a mathematical model 
may not have any physical meaning, so that we cannot observe them in the 
physical world. On the other hand, for a mathematical model, the (physical) 
distance between locations i and j may also not be well defined as in the 
physical world. 
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Fig. 3. Shape of the function p in Eq. (|2|). 

For the above reasons, in our experiments we chose to calculate dy based on 
the covariance matrix A itself. Suppose that A is a matrix with m rows (and 
the number of its rows is not less than the number of its columns), so that 



■l ' 



(3) 



where rj is the ith row of A. We define 



dij 



.T 



h/h 



(4) 



where l c is a length scale that is introduced to influence where the "cut-off" 
effect of the function p takes place 1 1. But note that here, l c does not have any 
physical meaning (i.e., it does not correspond to any quantity in the physical 
world). Instead, it is interpreted as a threshold of the statistical metric d^, 
which measures the (statistical) difference between the ith and jth elements 
of a random vector x. In [5J § 3.3.3.2], one of us showed through a numeri- 
cal example that covariance filtering conducted in this way can achieve the 



1 Note that we have assumed that the number of the rows of a matrix (not neces- 
sarily square) is not less than the number of its columns. If this is not the case, then 
it is suggested to choose the column vectors to calculate the distances dij in Eq. ([!]) 
instead. In this way, covariance filtering can be applied to non-square matrices like 
the cross covariance (when the dimension of the state space is not equal to that of 
the observation space). 



same effect as those reported in the literature by using the physical distances 
between different locations to construct the taper matrix. 

In [6], covariance filtering was conducted on the background covariance, the 
cross covariance, and the projection covariance (cf. Eq.(22) of j6]), since all of 
them are involved in updating the background to the analysis (cf. cf. Eq.(23) 
of [0]). No covariance filtering was conducted on the quantities other than the 
above three covariance matrices. 
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